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Abstract 

We derive expressions for the zeroth and the first three spectral moment sum rules for the 
retarded Green's function and for the zeroth and the first spectral moment sum rules for the 
retarded self-energy of the inhomogeneous Bose-Hubbard model in nonequilibrium, when the local 
on-site repulsion and the chemical potential are time-dependent, and in the presence of an external 
time-dependent electromagnetic field. We also evaluate these moments for equilibrium, where 
all time dependence and external fields vanish. Unlike similar sum rules for the Fermi-Hubbard 
model, in the Bose-Hubbard model case, the sum rules often depend on expectation values that 
cannot be determined simply from parameters in the Hamiltonian like the interaction strength 
and chemical potential, but require knowledge of equal time many-body expectation values from 
some other source. We show how one can approximately evaluate these expectation values for the 
Mott-insulating phase in a systematic strong-coupling expansion in powers of the hopping divided 
by the interaction. We compare the exact moments to moments of spectral functions calculated 
from a variety of different approximations and use them to benchmark their accuracy. 
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*Electronic address: jkf@physics.georgetown.edu 
^Electronic address: vturkows@mail.ucf.cdu 
+ Electronic address: hrkrish@gmail.com 



2 



I. INTRODUCTION 



The Bose-Hubbard model was originally introduced to describe the behavior of super- 
fluids in disordered environments, where superfluidity could be lost by phase fluctuations 
induced via disorder (Bose glass) and where potential-energy effects could localize the sys- 
tem into a Mott insulator pQ. While much work ensued on examining the properties of 
this model within the context of disordered superconductors, it was later proposed that this 
model could describe the behavior of bosonic atoms placed on optical lattices [2], which 
was experimentally realized soon thereafter [3J. Since the atoms in an optical lattice are 
also held to a small region of space via a trapping potential, the experimental systems are 
further complicated by being in an inhomogeneous environment. This is sometimes viewed 
as an advantage, as one can approximately map out the homogeneous phase diagram from 
the inhomogeneous system if the local-density approximation holds, but it also complicates 
much of the theoretical analysis, because the systems are no longer translationally invariant, 
so momentum is no longer a good quantum number. 

The model was extensively studied numerically even before the application to ultracold 
atoms was known. In this early work, the system was always assumed to be homogeneous. 
Monte Carlo simulations in one dimension [3] and two dimensions [5] , density matrix renor- 
malization group work [6], |7], as well as analytic approximations, like the strong-coupling 
approach [8], were carried out. After the experimental measurements of the Mott insulator 
were completed [3J, much further work ensued, culminating in the current state-of-the-art in 
quantum Monte Carlo simulations [9l-fTT]. in density-matrix renormalization group calcula- 
tions [32] and in experimental determinations of the phase diagram, including corrections to 
the local density approximation [T3HT6] . While much is known about the phase diagrams, 
less is known about the many-body spectral functions and self-energies of the Bose-Hubbard 
model. The many-body density of states (DOS) has been calculated via maximum entropy 
analytic continuation of quantum Monte Carlo (QMC) data [17J, via the variational cluster 
approach ( VCA) [18] , via the time-dependent density matrix renormalization group (DMRG) 
approach [19] and approximately with a bosonic version of dynamical mean-field theory in 
the strong-coupling limit [20] . 

Here, we describe how to derive exact sum rules for the spectral functions of the retarded 
Green's function (zeroth and first three moments) and of the retarded self-energy (zeroth 
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and first moment). We work with the most general case, which involves both spatially 
inhomogeneous and nonequilibrium (time-dependent) Hamiltonians, and also summarize the 
results for equilibrium and spatially homogeneous cases. This work complements previous 
work we have done on the nonequilibrium spectral moment sum rules for the Fermi-Hubbard 
model [221 - I24"] . which analyzed that model for the same general conditions. But unlike 
the Fermi-Hubbard model case, where all of the sum rules for retarded functions could be 
evaluated simply from parameters in the Hamiltonian, in the case for the Bose-Hubbard 
model, the spectral moment sum rules relate to many-body averages that are nontrivial, 
and need to be determined from some independent calculation in order to fully evaluate 
the different sum rules. This situation becomes more complicated for higher moments. 
But, these expectation values involve equal time expectation values for the nonequilibrium 
situation, and static expectation values for the equilibrium case, so they are simpler and 
more accurate to calculate than the spectral functions themselves. For example, we illustrate 
how to approximately determine them for the Mott-insulating phase with a strong-coupling 
perturbation-theory expansion. 

In Sec. II, we introduce the model and describe the techniques used to solve for the 
spectral moment sum rules and we develop the sum rules for the moments of the Green's 
functions and self-energies. We also discuss the equilibrium homogeneous case, and present 
the moments as functions of momentum. Numerical results comparing to different approx- 
imations for the spectral functions are performed in Sec. Ill along with a strong-coupling 
expansion for the expectation values needed to determine the values of the moments. Con- 
clusions follow in Sec. IV. 

II. FORMALISM FOR THE BOSE-HUBBARD MODEL 

The inhomogeneous nonequilibrium Bose-Hubbard model Hamiltonian can be written in 
the following general form in the Schrodinger representation: 

H{t) = -Y.tiAWh -Yl^A + \j2 U ^) b t b i( b l b i - !). (!) 

ij i i 

where fii(t) and Ui(t) can have arbitrary time dependence and can vary from site to site on 
the lattice. Similarly, the intersite hopping matrix —tij(t) can be time-dependent, such as 
in the case where there is an external electromagnetic field with arbitrary time-dependence 
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described by a vector potential A(r, t), and the bosons are charged, or in the case of cold 
atoms where the atoms move in a synthetic (Abelian) gauge field given by such a vector 
potential [21] • In this case, we would have a specific form for the time dependence of the 
hopping given by 



with q the "effective" charge of the bosons and R« the position vector for lattice site i. 

The general form for the Hamiltonian given in Eq. Q can be used to describe many 
different nonequilibrium and inhomogeneous situations. For example, it can describe bosons 
on an optical lattice that have the optical lattice potential modulated at some frequency, and 
thereby modifying the hopping, chemical potential and interaction as functions of time (as is 
done, for example, in modulation spectroscopy [25]); in this case, the local chemical potential 
would be described by the total chemical potential minus the local trapping potential at each 
lattice site. It can describe bosons on an optical lattice that have the trap either oscillated 
or impulsed to create dipole oscillations. It can also be used to describe a disordered Bose 
glass system (for a fixed distribution of local chemical potentials which represent the global 
chemical potential minus the diagonal on-site disorder). 

Similar to the fermionic case, one can define a generalized nonlocal spectral function for 
the boson retarded Green's function in the Heisenberg representation 



where we have introduced Wigner's time variables: the average time T = (ti + t2)/2 and the 
relative time t = t\ — t 2 (with t\ = T + t/2 the time at which the boson destruction operator 
is evaluated and t 2 = T — t/2 the time where the boson creation operator is evaluated). Note 
that we will often need to refer to the Green's function in both ways indicated in Eq. (|3]), 
in terms of the times the operators are evaluated t\ and t 2 , or in terms of the average 
and relative times. While we do not introduce a new notation for the reexpressed Green's 
function, the context for whether the time variables are the variables where the operators 
are evaluated or are the average and relative times is always obvious by the context of the 







in the following way: 




(4) 
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equations and whether the time variables have a subscript or not. The square brackets are 
used to denote a commutator. The angular brackets denote a statistical average with respect 
to the initial equilibrium distribution of the system in the distant past, with a density matrix 
given by p = exp[—/3H(t — » —oo)]/Z and Z = Trexp[— (3H(t — > — oo)]. Note that T denotes 
the average time, not the temperature, which is denoted by 1//3. The time dependence of 
the operators is with respect to the full Hamiltonian and is expressed in the Heisenberg 
representation via an evolution operator. 

We further define the nth spectral moment of the corresponding spectral function in 
Eq. Q to be 

POO 

p%{T) = / duu n AfAT,u). (5) 



This spectral moment is connected to the bosonic Green's function in the following way: 

1 



Tm 



7T 



dco 



an 



dt n lJ 



From this equation, one can obtain the following formula: 

& 



Tm 



(6) 



(7) 



J t=o+ 



dt n 13 

It is possible to show that the terms proportional to the time derivatives of the theta function 
do not contribute to the moments. Then, by using the Heisenberg equation of motion for 
the operators, one can obtain the following expressions for the zeroth and the first three 
spectral moments of the Green's function: 



(MT),b](T)}), 



([L\{T)MT)\) - ([kinL^T)}) 



+ 



([L\(T),b](T)]) - 2([L\(T),L'b](T)]) + ([6,(T), L%](T)}) 

; ,(nH'(T)],b](T)]) + [b](T),H'(T)]]}_ 

([L\(T), b](T)\) - S([L%(T), L x b)(T)\) + m\{n L 2 b](T)}} 
(MT),L 3 b](T)}}_ 

q mHnH\T)],H(T)],b](T)]) - ([[[b t (T),H(T)},H'(T)},b](T)}} 
-3([[[h(n H'(T)],b](T)], H(T)}} + 3([[\b t (n H ( T )l b](T)] 7 H'(T)]) 
-(UT) ,[[b](T),H'(T)],H(T)]) ~ 2([b t (T),[[b](T),H(T)},H'(T)}) 
i(nH"(T)],b](T)]) - ([bi(T), [bj(T), H"(T)]])~ 



(8) 
(9) 

(10) 



% 

+ - 



(11) 



where L n = \..\[0,H],H]...H] is the operation of commutation performed n times, and 
H'(T) = dH(T)/dT and H"(T) = d 2 H{T)/dT 2 are the first and the second derivative of 
the Hamiltonian in the Schrodinger picture with respect to time which is then evaluated in 
the Heisenberg picture at the average time T (i. e., the derivative is taken before going to 
the Heisenberg picture). In equilibrium these derivatives are zero, but in the nonequilibrium 
case they can be finite due to the explicit time-dependence of the Hamiltonian parameters. 
Evaluating the commutators is a straightforward, but tedious exercise. Doing so results in 
the following expressions: 

t4j(T) = 6 ijt (12) 
/if^T) = —tij(T) + 2U t (T) ni (T)5 l3 , (13) 

i 

~Ui{T)rii {T)5ij + 3U?(T)(n 2 (T))5ij, (14) 



^3ij( T ) - -y]tjltlmtmj(T) 

l,m 

+2U ini J2 Uitij{T) + 2 J2 tiiUmiiiAT) + 2 ^ /;,/; ; r / // ; (V) 
i i i 

+2Ui8ij (tniin(b\b n ) + tut^blk) - 2i u tni{bibi)) (T) 

l,n 

+U 2 n^(T) + tijU 2 rij(T) - 30?<n?)f o -(T) - 3tyE/J >(n)){T) 
-AUiUjUjirnnjXT) - UjijiUiib^bi^+iiiUfrH^Sij 



+2UA,J2 Ul [tiMb^+tiMk)) (Tj-AUfSij^tiMbiKT) 
i i 

-2U t 5 l3 J2 Ul (tuiblbmij+tuiniblbi)^ (T) + 6U 2 S io *«< 6 I W( T ) 



i 

-?>5 l3 Ut(n 2 )(T) + SijUfrii (T) + A5 l3 Uf (n\)(T) 
+ 1 E - tiAj) ( T ) - IWijUj ~ U 3 U' 3 )n 3 + {Ufa - Uj tJ )n t ](T) 

where m(T) = (b\(T)bi(T)), -Uj(T) = -Uj(T) - <%//i(T), and the symbol (T) reminds 
us that both the parameters in the Hamiltonian and the operator expectation values are to 
be evaluated at the average time T in the expressions for the moments. Note that we have 
expressed the operator expectation values in the most compact form rather than subtracting 



off the average values to show the effects of correlations about the average values. It is a 
simple, but tedious, exercise to convert to other forms if desired. 

The expressions for the retarded self-energy moments can be obtained by using the Dyson 
equation, which connects the retarded Green's function and self-energy 

Ggfa.ta) = Gf\t u t 2 ) + E J dt * f dUGf l { "\t 1 M)^UhM)Gi 1 {UM), (16) 



where Gff°\t 1 , t 2 ) is the noninteracting Green's function (in the case U = 0), and S^(t 1) i 2 ) 



is the retarded self-energy. The first step is to rewrite Eq. (16) in a combined frequency- 
average time representation: 

Gg(7» = G™(T,cu) + J2 f dT I di ! dn ! dve-^ l e lvf 

lm J J 

x G^^T+| + |,w + n + 0S^(T + f |W + 2n)G^^T+|-^ W + n-0. 

(17) 

where the internal time variables with the overbars should not be confused with the effective 
hopping matrix used to describe the Green's function moments. 

One can represent the Green's function and the self-energy at large frequency in terms 
of an asymptotic l/oo expansion: 

OO ~ R (rp\ 

n=0 

= E i^r- + ^(T," = oo). (20) 

n=0 

where fi^ n (T) denotes the corresponding Green's function moments for the noninteracting 
Green's function (which has U — 0) and C^AT) are the moments of the retarded self-energy, 



defined via 



1 f°° 

C^{T) = --\mj duu n ^(T,co). (21) 
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The large-u; expansions in Eqs. (18)-(20) can be obtained by using the following spectral 



identities (valid for all retarded functions that decay rapidly enough for large relative time): 

^^-'f.M, (22) 

J IX J_ 00 OJ - Oj' + l0 + 

At large u> the Green's function and self-energy must be real, or the spectral moments would 
diverge. The self-energy expansion in the last equation contains a frequency-independent 
term £^(T, u> = oo), which corresponds to the constant (Hartree-Fock) term in the self- 
energy. 

Then, one inserts these expansions into Eq. (17) and considers separately all terms which 
have the same order in (1/u)). Hence, it is necessary to expand all the functions under the 
integrals in powers of (1/w). This leads the following equations which connect the Green's 
functions and self-energy spectral moments: 

/< = = hi, (24) 
Vuj = fiiij + Y VoiAP, w = oo)/i^., (25) 



Lm 



R ~ R ^ ^ ~ R v 1 R St* \ R ^ ^ ~ R s~iR R 

Lm l,m 

+ fiiA(T, uj = oo)n* mj , (26) 



Lm 



,.R — r,R i\^r, R ^p R (T,, — ^\,, R _i_ V* r, R n R ,, R _u V* r, R n R ,, R 

rZij ~ nij "T" / .H'Oil^lmK 1 ^ ~ °°JH'2mj ' 2-U r'Oil^QlmH'lmj ' / 4 h^Oil^ llmH'Omj 
l,m Lm l,m 

+ Y fiiiPL(T, uj = oo)/if mj + Y, VuiCLvLj + Y ~^L(T, u = oo)/i K mj .(27) 

l,m l,m l,m 

We have suppressed the average time dependence of all of the spectral moments to save 
space. 

From these equations, one can obtain the results for the retarded self-energy moments 
by using the results for the retarded Green's function moments derived above for both finite 
and vanishing U: 

Sg(T, oj = oo) = 2U l n l 5 ij iT), (28) 
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C * (T) = -UfnAAT) ~ Wtn%i{T) + 3U?5 tl (n*)(T), (29) 
= 2U i 5 ij J2(iiiiin(b\bn)+t H ini(blb i ) - 2i il inM l b l )) (T) 

l,n 

■ U'J ,jl' jdi.ii, - (ninj))(T) 

+<S«tf?(4n? + 8n? - 12^) - 3(n?> + n, + 4<n J 3 ))(T) 
-U^Ui^hb^iT) + UiUfmSijlT) 

+2UAj U i (tn(blbi) + t^&I&i)) (T) - 4C/^.. J^tuiblk)^) 
i i 

-2U t 5 tJ J2 U i (tuiblbm^+tuimblb^ {T) + ml5 t] Y J Ub%n l )(T) 
i i 

-\mP 3 - UjUfc + (Ufa - mjhiKT) - i^P n ,(T)%, (30) 

where, once again, the (T) symbol is to remind us of the average time dependence of both 
the parameters of the Hamiltonian and of the operator expectation values. 

The expressions for the spectral moment sum rules are complicated in general. We want 
to summarize these results for the case of a homogeneous system in equilibrium, where 
there is no time dependence to the Hamiltonian, and where the hopping is the same value 
between all nearest neighbors, the chemical potential is uniform, as is the interaction energy 
U . In this case, we can express the moments in momentum space, with respect to the 
time-translation-invariant momentum dependent Green's function 

G«(t) = -ie(t)([b*(t),bt(0)]), (31) 

and self-energy 

S k » = G™H- 1 -G^H- 1 , (32) 

where the symbol u> is used for the Fourier transform from time to frequency space, which 
is used because the Hamiltonian is time independent in equilibrium. The operator b^ is the 
momentum-space destruction operator, which satisfies &k = X)i h exp[— ik • R^/v^/V, with a 
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corresponding formula for We find the Green's function moments become 




1 



(33) 




£k — A* + 2Un, 



(34) 




= (e v -n + 2Un) 2 + U 2 (3(n 2 ) - An 2 - n) , 

= (e k - n + 2f/n) 3 + f/ 2 Zt ii+5 (^6^ l+5 + 6t +<5 6 ini ) 

+ f/ 2 e k (6(n 2 ) - 12n 2 - 2n + 4(n^ i+5 ) + - 3£/ 2 /i (3(n 2 ) - 4n 2 

+ £/ 3 (4(n 3 ) - 8n 3 - 3(n 2 ) + n) . 



(35) 




(36) 



Here, the symbol 5 denotes the nearest neighbor translation, (nin i+ s) denotes the nearest- 
neighbor static density-density expectation value and —t ii+ s is the nearest-neighbor hopping 
matrix element. Since the system is homogeneous, the expectation values and the hopping 
matrix element are independent of which nearest neighbor translation 5 is chosen. The high- 
frequency constant for the self-energy is E(tu — > oo) = 2Un, and the self-energy moments 
become 



The zeroth and first Green's function moments and the constant of the self-energy do 
not require any expectation values besides the filling. The second Green's function and 
zeroth self-energy moments require just one correlation function (n 2 ) — n 2 , while the higher 
moments require increasingly more (and more complex) correlation functions. 

III. NUMERICAL RESULTS 

There are not too many techniques which can accurately determine the Green's function 
and self-energy of the Bose Hubbard model for real frequencies. To date, the main non- 
perturbative methods that have been tried include quantum Monte Carlo plus a maximum 
entropy analytic continuation |17j . the variational cluster approach [T5], the time- dependent 
density matrix renormalization group [19], and a strong-coupling version of bosonic dynam- 
ical mean-field theory [20] . The dynamical mean- field theory calculation is most accurate in 




U 2 (3(n 2 ) - An 2 - n) , (37) 
U 2 Zt vl+s ( ni b\b i+s + fet^rii) + t/ 2 e k (4(n,n 4+5 ) + (bf +5 b 2 ) - 4n 2 ) 
U 2 fi (S(n 2 ) + An 2 + n) + U 3 (4(n 3 ) + 8n 3 - 3(n 2 ) - 12(n 2 )n + An 2 + n) . (38) 



a 



iR 
Ik 



+ 
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large dimensions, and we will not consider it further here. Both the QMC and VCA can be 
performed in any physical dimension, the density matrix renormalization group is limited to 
one dimension, which is where we will focus our attention first. The quantum Monte Carlo 
approach can also calculate different static expectation values, so using this data will allow 
us to completely determine the different moments. 

But, in general, it would be nice to have alternative methods to at least approximate 
the value of the spectral moment sum rules. In the case of a Mott insulator with a small 
hopping, one can use strong-coupling perturbation theory in the hopping to calculate the 
different expectation values. The Mott phase for a homogeneous system in equilibrium with 
vanishing hopping is given by 

A ("0' 



n 70 



IIVrl >' (39) 



i=i Vn< 



where |0) is the vacuum state and the Mott phase has an average density of n. The energy 
of this state is E® = [—fin + Un(n — 1)/2]JV. Since this state is nondegenerate, standard 
nondegenerate Rayleigh-Schrodinger perturbation theory can be used to find the wavefunc- 
tion. Proceeding in the canonical fashion, where the perturbed state (n\ satisfies (n\n) = 1, 
we find that the perturbed wavefunction satisfies 

\n) = |n) + ® n „ V\n) + ^\ V ^\ V\n) , (40) 
E° n -H E°-H E® — Hq 

through second order in the perturbation V, since the first-order shift in the energy 
o(n\V\n) = vanishes (the perturbation V is the hopping term in the Hamiltonian). Here 
Q n = I — \n)o o(n\ is the projector onto the space perpendicular to the unperturbed wave- 
function |n)o- A straightforward calculation of the overlap of the perturbed wavefunction 
is 

Zt 2 ( t 4 \ 

(n\n) = 1 + — n(n + l)N + ^— J . (41) 

Here Z is the number of nearest neighbors and t is the magnitude of the nearest-neighbor 
hopping matrix element (we assume we are on a bipartite lattice so all odd order terms in 
the perturbation V vanish). 

It is now a straightforward exercise to find the expectation values needed for the different 
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TABLE I: Comparison of expectation values calculated in strong coupling perturbation theory 
versus quantum Monte Carlo simulation. The case considered is a one-dimensional Mott insulator 
with n = 1, t = 0.05U, and low temperature. 



Expectation value 


strong coupling 


quantum Monte Carlo 


(n 2 ) 


1.02 


1.02 ±0.0004 


<*?> 


1.06 


1.06 ±0.001 


{nim+d) 


0.99 


0.9902 ± 0.0002 




0.20 


0.2022 ± 0.0004 




0.01 


not calculated 



moments. We have 

{n\rii\n 



(n\n) 
(n\n 2 \n 



{n\n) 
{n\n 3 \n 



(n\n) 
(n\nin i+s \n 



(n\n) 



(n\n t b\b i+5 \n 



{n\n) 



\n 



n, 



2 2Zt 2 . . „{t 4 \ 



n 



QZt 2 
2t 



n 2 (n + 1) + O 



n 2 - j^nin + 1) + O 



(n\n) 



2t 
U 

u 2 



u 2 

n 2 (n + 1) +0 
n(n + l)[n(n - 



u 3 

n- 



\{n-l){n + 2)]+oU Ji 



(42) 
(43) 
(44) 
(45) 
(46) 
(47) 



Here 5 denotes a nearest neighbor translation, so that i + 5 is a nearest neighbor site of site 
% (the expectation values are independent of which neighbor because the problem we are 
solving is homogeneous). 

We compare the accuracy of the above expectation values with ones calculated directly 
via a quantum Monte Carlo approach in Table |T] for a Mott insulator on a one-dimensional 
lattice with t/U = 0.05, n — 1, and low temperature. One can see that the accuracy is 
excellent for the strong coupling perturbation theory in this parameter regime. 

Using the strong-coupling perturbation theories, we find the momentum-dependent re- 
tarded Green's function moments approximately become 

6Zt 2 ' 



^ 2k 



(e k - 11 + 2Un) 2 + U 2 n(n + 1) 
13 



1 



U 2 



(48) 



^ 3k 



e k - n + 2Un) 6 + AUZt 2 n 2 (n + 1) - 2U 2 e k n(n + 1) 



x 



1 - 



2Zt 2 



t 2 



u 2 u 2 



[n(n + l) - -(n- l)(n + 2)] 



+ 3£/ 2 /m(™ + 1) 1 



IT 2 " 



£/ 3 n(n + l)(4n 



It 2 " 



(49) 



to order t 4 /[/ 4 . The self-energy moments approximately become 

6Zt 2 ' 



°0k 



-U z n (n + 1 1 - 



£/ 2 



(50) 



and 



= 4UZt 2 n 2 (n + 1) - t 2 e k n(n + 1)[8 - n(n + 1) - -(n - l)(n + 2)] 



"tT 2 " 



+ [/ 3 ri(ri + 1)1 



1 

2 

6Zt 2 

IT 2 " 



(51) 



also to order t A /U 4 . 

Now we are ready to compare the accuracy of different exact methods in calculating the 
many-particle density of states for the Bose-Hubbard model. Our first test case is the Mott 
insulating phase in the one-dimensional model with t/U = 0.05, n = 1, and fi/U = 0.5 
In Fig. [TJ we plot the local density of states for the three different methods that have 
been used for this problem: (i) the VCA at zero temperature with a Lorentzian broadening 
of r) = 0.03 and rj = 0.002 [18]; (ii) QMC simulation plus maximum entropy analytic 
continuation [T7j, where the calculations are performed at a temperature (3 = 192; and (iii) 
time- dependent density matrix renormalization group calculations [19] at zero temperature 
with open boundary conditions and a Lorentzian broadening of rj = 0.04. One can see that 
there is a significant discrepancy between these different curves (in a pointwise fashion) but 
as we will see the moments are all quite close to one another. This tells us that the density 
of states is quite sensitive to the broadening chosen, and it is difficult to tell which of these 
different techniques is most accurate (although, the quantum Monte Carlo technique uses 
the most unbiased algorithm to determine the density of states). It is also apparent that 
one properly sees the correct gap in the density of states only with the methods that use 
the least broadening, as expected. 

Now we move on to the task of comparing the different spectral moment sum rules. We 
begin with the zeroth moment sum rule of the retarded Green's function, which essentially 
tests whether the system has conserved the correct number of states in the given calculation. 
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FIG. 1: (Color online.) Local density of states for the one-dimensional Bose Hubbard model 
in the Mott-insulating phase with n = 1. The parameters are t/U = 0.05 and \xj\J = 0.5. The 
energies are measured in units of U. We compare the variational cluster approach with two different 
broadenings (red and blue) to the quantum Monte Carlo plus maximum entropy approach (purple) 
and to the density matrix renormalization group approach (orange). The inset zooms in on the 
region just above 2U, where the variational cluster approach has some structure which is needed 
to get high precision to the different moment sum rules. The data for the other two methods is 
cutoff before this frequency. 

For the variational cluster approach, since it determines the Green's function as a set of 
delta functions and weights, we perform the integration of the moments via summing the 
relevant weights of the delta functions rather than introducing any artificial broadening 
into the calculation. Doing this, appears to greatly improve the accuracy of the spectral 
moments themselves. In Fig. [2] (a) we plot the exact result against the three different 
approximation techniques. On the graph, one cannot see any error between the VCA and 
the exact result. The quantum Monte Carlo is about a half a percent error, while the density 
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FIG. 2: (Color online.) Spectral moments of the retarded Green's function as a function of momen- 
tum for the one-dimensional Bose Hubbard model in the Mott-insulating phase with n = 1. The 
parameters are t/U = 0.05 and fi/U = 0.5. We compare the exact result (black) to the variational 
cluster approach (red dashed line), to the quantum Monte Carlo plus maximum entropy approach 
(purple), and to the density matrix renormalization group approach (orange). Panel (a) is the 
zeroth, panel (b) the first, panel (c) the second and panel (d) the third moment. The accuracy of 
the VCA is so good, one cannot see any deviation from the exact result in panels (a) and (b). The 
quantum Monte Carlo and density matrix renormalization group results have higher errors. 

matrix renormalization group results are about three and a half percent error. 

The first retarded Green's function moment is plotted in Fig. [2] (b). Here the VCA and 
the QMC have small errors, with the former being less than 0.1% and the latter on the order 
of a few percent. The density matrix renormalization group result has the right shape, but 
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appears to be systematically shifted off of the correct result causing about a 7% error. 

The second retarded Green's function moment is plotted in Fig. [2] (c). This is the first 
moment that depends on a correlation function. Once again, the VCA has superior accuracy, 
and the density matrix renormalization group results look like they are systematically shifted 
from the correct answer (so much so that at small momentum they have the wrong sign for 
the moment). 

Finally, in Fig. [2] (d), we plot the third moment sum rules. Here the deviations of all 
of the approximations are larger, but the VCA is clearly superior. One might ask why the 
VCA appears to be so much better than the other two techniques, at least when we compare 
the moment sum rules? We believe the answer to this lies in the inset in Fig. [I] There one 
can see that the VCA has some spectral weight at high energy, corresponding to higher 
on-site occupation of bosons. The QMC and density matrix renormalization group results 
are cut off at lower frequencies, so they do not have this extra feature. This feature becomes 
increasingly important for higher moments, since the integrands are weighted more and also 
for even moments, since it can modify the cancellation that occurs between the positive and 
negative branches of the density of states. In this sense, the moment sum rules can be a 
very delicate test of the accuracy of the different numerical calculations. In addition, the 
VCA, being based on a strong-coupling approach, is more accurate for large U . We would 
expect other techniques to become more accurate as we move closer to the critical point and 
beyond. 

We also want to test the spectral moment sum rules of the self-energy. Here we have to 
now further process the data, as one cannot compute the self-energy solely from the density 
of states or the momentum-dependent spectral functions. Instead, we use a Kramers-Kronig 
relation on the spectral functions to determine the momentum-dependent retarded Green's 
functions (the imaginary part is just — 7r times the spectral function). Then we use Dyson's 
equation to extract the self-energy. As a test, we compare the constant term of the real 
part of the self-energy to its exact result. For the quantum Monte Carlo and for the density 
matrix renormalization group data, the extent of our data is too small in frequency for us to 
properly reach the limit where we can accurately extract the constant, but the error in the 
constant is less than 15%. Once we have the imaginary part of the self-energy, we simply 
integrate it times the appropriate power of frequency to see the sum rule. For the VCA, we 
can no longer do this with the delta function representation, so we instead use the smaller 
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broadening r] = 0.002£7 data, which then leads to some noisy fluctuations in the integrated 
self-energy moments. But the total noise level is not too high. 
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FIG. 3: (Color online.) Zeroth (a) and first (b) spectral moments of the retarded self-energy as a 
function of momentum for the one- dimensional Bose Hubbard model in the Mott-insulating phase 
with n = 1. The parameters are t/U = 0.05 and ji/U = 0.5, with energies measured in units of 
U. We compare the exact result (black) to the variational cluster approach (red dashed line), to 
the quantum Monte Carlo plus maximum entropy approach (purple), and to the density matrix 
renormalization group approach (orange). Deviations are visible for all approximations. 

We plot the zeroth spectral moment sum rule for the retarded self-energy as a function 
of momentum in Fig. [3] (a). The sum rule itself is independent of momentum, even though 
the self-energies do vary with momentum. The VCA has errors at the 0.5% level. The 
density matrix renormalization group has errors about ten times larger, and once again, 
there appears to be a systematic error in that data pushing it to slightly lower values. 

Finally, we plot our last one-dimensional result, the first moment spectral sum rule of the 
retarded self-energy versus momentum in Fig. [3] (b). Here we see similar results as with the 
zeroth moment, perhaps with somewhat larger errors. 

The conclusion of this work on the one-dimensional example is that the spectral moment 
sum rules for the retarded Green's functions and for their self-energies provides useful data 
to help us predict the accuracy of different numerical calculations. While they cannot 
provide us with sufficient data to determine what the appropriate widths of different spectral 



18 



features should be in the Green's functions by examining pointwise values of the spectral 
functions, they do tell us important information about the weight under the curves and 
of their respective shapes. Indeed, the higher moments, are particularly sensitive to small 
weight structures at high energy and could help identify whether approximations are cut off 
at too small a frequency and missed some higher energy spectral weight. 

For two dimensions, we only have data from the VCA and from the QMC. We have ana- 
lyzed all of the moments in a similar fashion to what was done for the one-dimensional case, 
and we summarize our results in a series of figures. We choose parameters corresponding to 
the second Mott lobe with n — 2, and with a hopping that lies inside the lobe but close to 
the tip (t/U = 0.03 and fi = 1.45). In this case, we might expect to see larger deviations 
from the sum rules. The quantum Monte Carlo results are run at a low temperature (3 = 80. 
In Figure |4j we plot the exact results for the moment sum rules (as evaluated from the 
strong-coupling perturbation theory described above) versus the VCA and quantum Monte 
Carlo results from a maximum entropy analytic continuation. Since the momenta are now 
distributed through the two-dimensional Brillouin zone, we show plots for momenta along 
the high-symmetry lines of the triangle that runs from the origin at the T point (r = 0) 
to the M point along the diagonal direction [M = (tt,tt)], to the X point along the axial 
direction [X = (jr, 0)]. One can see for the low moments, the VCA approximation continues 
to work extremely well (once again, the VCA is evaluated as sums over the delta functions, 
and uses no broadening, which is why the moments are so accurate). But even in this case, 
we do start to see deviations for higher moments, and they are larger than they were for the 
one-dimensional case; this is expected due to the proximity to the tip of the Mott lobe at 
t = 0.035 (26]. The QMC results, on the other hand, do show the right trends, but have a 
much larger quantitative error to the moments. 

We continue with plots of the zeroth and first moments of the self-energy in Fig. [5] Here, 
we must go through the inversion procedure described above to extract the self-energy from 
the data for the Green's function. Hence, for the VCA, we must broaden the delta functions 
to create a smooth functional form for the Green's function. We do this with a narrow 
broadening parameter to try to preserve the accuracy. This gives rise to the amplitude of 
the oscillations in the moment data due to oscillations in the Green's function and then in 
the self-energy. Once again, we do see the right trends in the data, but here the accuracy 
is fairly poor, especially for the QMC data. The VCA data tends to have the right average 
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FIG. 4: (Color online.) Spectral moments of the retarded Green's function as a function of mo- 
mentum for the two-dimensional Bose Hubbard model on a square lattice in the Mott-insulating 
phase with n = 2. The parameters are t/U = 0.03 and fi/U = 1.45. We compare the exact result 
(black) to the variational cluster approach (red dashed line), and to the quantum Monte Carlo 
plus maximum entropy approach (purple). Panel (a) is the zeroth, panel (b) the first, panel (c) 
the second and panel (d) the third moment. The accuracy of the VCA is so good, one cannot see 
any deviation from the exact result in panels (a) and (b). The quantum Monte Carlo results have 
higher errors. The symmetry lines over which the data is plotted are shown with the schematic 
triangle. 
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FIG. 5: (Color online.) Zeroth (a) and first (b) spectral moments of the retarded self-energy as 
a function of momentum for the two-dimensional Bose Hubbard model on a square lattice in the 
Mott-insulating phase with n = 2. The parameters are t/U = 0.03 and n/U = 1.45, with energies 
measured in units of U. We compare the exact result (black) to the variational cluster approach 
(red line), and to the quantum Monte Carlo plus maximum entropy approach (purple). Deviations 
are visible for all approximations indicating the self-energy is not determined as accurately here. 

behavior, but it's quantitative average is incorrect. By looking at the traces of the self-energy 
itself, there are likely two reasons for the discrepancy. The first, is that the frequency cutoff 
might be too low, and the second is that the self-energy seems to have regions of frequency 
where it has strong frequency dependence, and this might not be properly captured by the 
approximations. 

IV. CONCLUSIONS 

In this work, we have derived exact formulas for the first few spectral moments of the 
Bose-Hubbard model through third order for the retarded Green's function and through 
first order for the retarded self-energy. The results we derive are quite general, holding for 
inhomogeneous systems and for systems that have time dependence to the parameters in the 
Hamiltonian. Sum rules can be quite useful in benchmarking different approximations, be- 
cause their results are exact. One challenge with the work here, is that for the bosonic case, 
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the moments depend on correlation functions that must be determined for the interacting 
system, unlike in the fermionic version, where many of the correlation functions become triv- 
ial to evaluate. But, in the limit of strong coupling, for the Mott phase, these moments can 
be systematically found in a strong-coupling expansion, which appears to be quite accurate 
when compared to QMC results. We concluded this work with numerical calculations for 
translationally invariant systems in equilibrium, where we could benchmark the accuracy of 
different numerical results. Because we did this in the strong-coupling region, it comes as 
no surprise that the VCA turned out to be the most accurate approach, indicating that it is 
faithfully producing the moments of the spectral functions. It is much more difficult for us 
to determine the point-wise accuracy of the spectral functions, though. We did see that the 
numerical calculations appear to work best in one-dimension, where the moment sum rules 
are most accurate. 

We hope that this work will be used by others for benchmarking purposes of numerics 
and possibly for understanding qualitative changes in spectra as parameters change due to 
the constraints given by the sum rules. As nonequilibrium techniques are developed for 
interacting bosonic systems, it will also be interesting to use these results for benchmarking 
of those calculations, since exact results for nonequilibrium systems are quite rare. 
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